! include 'mkl_pardiso.f90'

module MySparseOperator
	implicit none
contains

	subroutine mso_zgmres_precondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres, prefun)
		!*******************************************************************************
		!  Content:
		!  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
		!*******************************************************************************
		implicit none
		interface
			subroutine prefun(x,u)
				complex(kind=8), dimension(:), intent(in) :: x
				complex(kind=8), dimension(:), intent(out) :: u
			end subroutine prefun
		end interface
		integer, intent(in) :: maxiter
		real(kind=8), intent(in) :: maxrelres
		complex(kind=8), dimension(:), intent(in) :: a
		integer, dimension(:), intent(in) :: ia, ja
		complex(kind=8), dimension(:), intent(in) :: rhs
		complex(kind=8), dimension(:), intent(out) :: computed_solution
		integer, intent(out) :: itercount
		real(kind=8), intent(out) :: relres
		

    	integer :: n, zn
    	integer :: ipar(128)
    	real(kind=8) :: dpar(128)
    	real(kind=8), dimension(:), allocatable :: tmp, b, drhs, dcs
    	complex(kind=8), dimension(:), allocatable :: residual
    	
    	integer :: rci_request, i
    	real(kind=8) :: dvar, rhsm2
    	real(kind=8), external :: dznrm2
    	integer :: err
    	
    	err = 0

    	zn = size(ia) - 1
    	n = 2*zn	! the size of the problem
    	rhsm2 = dznrm2(zn,rhs,1)

    	allocate( b(n))
    	allocate( residual(zn))
    	allocate( drhs(n))
    	allocate( dcs(n))

		!---------------------------------------------------------------------------
		! Initialize variables and the right hand side through matrix-vector product
		!---------------------------------------------------------------------------
    	call dcopy(zn, dreal(rhs), 1, drhs, 1)
    	call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

    	dcs = 1.d0

    	ipar(15) = min(1500,maxiter)
    	allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

    	call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999
    	
    	ipar(8) = 1
    	ipar(15) = min(1500,maxiter)
    	ipar(5) = ipar(15)
    	ipar(11) = 1
    	dpar(1) = maxrelres

    	call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999

    	1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	
    	if (rci_request .eq.0) go to 3

    	if (rci_request .eq. 1) then
    		computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),	&
    				&	tmp(ipar(22)+zn:ipar(23)-1) )
!     		call zcopy(zn,dcmplx(tmp(ipar(22):ipar(22)+zn-1),tmp(ipar(22)+zn:ipar(23)-1)),1,computed_solution,1)
    		call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
    		call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
    		call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)		
            go to 1
    	end if

    	if (rci_request .eq. 2) then
 			! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
			! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
			! compute the current true residual via mkl (sparse) blas routines
			computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
        	computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),	&
        			&	tmp(ipar(22)+zn:ipar(23)-1) )
        	call prefun( computed_solution, residual )
        	call dcopy(zn, dreal(residual), 1, tmp(ipar(23)), 1)
        	call dcopy(zn, dimag(residual), 1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

    	3 ipar(13) = 0
    	call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
    	computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2
    	
    	if (allocated(dcs)) deallocate(dcs, stat=err)
    	if (err /= 0) print *, "dcs: Deallocation request denied"
    	if (allocated(drhs)) deallocate(drhs, stat=err)
    	if (err /= 0) print *, "drhs: Deallocation request denied"
    	if (allocated(tmp)) deallocate(tmp, stat=err)
    	if (err /= 0) print *, "tmp: Deallocation request denied"
    	if (allocated(b)) deallocate(b, stat=err)
    	if (err /= 0) print *, "b: Deallocation request denied"
    	if (allocated(residual)) deallocate(residual, stat=err)
    	if (err /= 0) print *, "residual: Deallocation request denied"
    	return
    	

    	999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
    	call mkl_free_buffers
    	stop 1
	end subroutine mso_zgmres_precondition




	subroutine mso_zgmres_nonprecondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres)
		!*******************************************************************************
		!  Content:
		!  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
		!*******************************************************************************
		implicit none
		integer, intent(in) :: maxiter
		real(kind=8), intent(in) :: maxrelres
		complex(kind=8), dimension(:), intent(in) :: a
		integer, dimension(:), intent(in) :: ia, ja
		complex(kind=8), dimension(:), intent(in) :: rhs
		complex(kind=8), dimension(:), intent(out) :: computed_solution
		integer, intent(out) :: itercount
		real(kind=8), intent(out) :: relres
		

    	integer :: n, zn
    	integer :: ipar(128)
    	real(kind=8) :: dpar(128)
    	real(kind=8), dimension(:), allocatable :: tmp, b, drhs, dcs
    	complex(kind=8), dimension(:), allocatable :: residual
    	
    	integer :: rci_request, i
    	real(kind=8) :: dvar, rhsm2
    	real(kind=8), external :: dznrm2
    	integer :: err
    	
    	err = 0

    	zn = size(ia) - 1
    	n = 2*zn	! the size of the problem
    	rhsm2 = dznrm2(zn,rhs,1)

    	allocate( b(n))
    	allocate( residual(zn))
    	allocate( drhs(n))
    	allocate( dcs(n))

		!---------------------------------------------------------------------------
		! Initialize variables and the right hand side through matrix-vector product
		!---------------------------------------------------------------------------
    	call dcopy(zn, dreal(rhs), 1, drhs, 1)
    	call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

    	dcs = 1.d0

    	ipar(15) = min(1500,maxiter)
    	allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

    	call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999
    	
    	ipar(8) = 1
    	ipar(15) = min(1500,maxiter)
    	ipar(5) = ipar(15)
    	ipar(11) = 0
    	dpar(1) = maxrelres

    	call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999

    	1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)
    	
    	if (rci_request .eq.0) go to 3

    	if (rci_request .eq. 1) then
    		computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),	&
    				&	tmp(ipar(22)+zn:ipar(23)-1) )
    		call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
    		call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
    		call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)		
            go to 1
    	end if

    	if (rci_request .eq. 2) then
 			! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
			! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
			! compute the current true residual via mkl (sparse) blas routines
			computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
        	call dcopy(n,tmp(ipar(22)),1, tmp(ipar(23)), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

    	3 ipar(13) = 0
    	call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
    	computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2
    	
    	if (allocated(dcs)) deallocate(dcs, stat=err)
    	if (err /= 0) print *, "dcs: Deallocation request denied"
    	if (allocated(drhs)) deallocate(drhs, stat=err)
    	if (err /= 0) print *, "drhs: Deallocation request denied"
    	if (allocated(tmp)) deallocate(tmp, stat=err)
    	if (err /= 0) print *, "tmp: Deallocation request denied"
    	if (allocated(b)) deallocate(b, stat=err)
    	if (err /= 0) print *, "b: Deallocation request denied"
    	if (allocated(residual)) deallocate(residual, stat=err)
    	if (err /= 0) print *, "residual: Deallocation request denied"
    	return
    	

    	999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
    	call mkl_free_buffers
    	stop 1
	end subroutine mso_zgmres_nonprecondition



	subroutine mso_dcoo2csr(rowind,colind,acoo,acsr,ia,ja,m,nnz)
		integer, intent(in) :: m,nnz
		integer, intent(in) :: rowind(nnz), colind(nnz)
		real(kind=8), intent(in) :: acoo(nnz)
		real(kind=8), intent(out) :: acsr(nnz)
		integer, intent(out) :: ia(m+1), ja(nnz)

		integer :: job(8), info

		job(1) = 2	! converted to CSR
		job(2) = 1	! one-based indexing for CSR
		job(3) = 1	! one-based indexing for COO

		job(6) = 1
		call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
		job(6) = 2	
		call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
	end subroutine mso_dcoo2csr

	subroutine mso_dcsradd(trans,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic)
		character, intent(in) :: trans
		integer, intent(in) :: m,n, ia(*), ja(*), ib(*), jb(*)
		real(kind=8), intent(in) :: beta, acsr(*), bcsr(*)
		real(kind=8), dimension(:), intent(out) :: ccsr
		integer, dimension(:), intent(out) :: jc, ic

		integer :: info, nzmax

		call mkl_dcsradd(trans,1,3,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic,nzmax, info)
		call mkl_dcsradd(trans,2,3,m,n, acsr,ja,ia, beta, bcsr,jb,ib, ccsr,jc,ic,nzmax, info)
	end subroutine mso_dcsradd

	subroutine mso_dcsrdltRC(a,ia,ja,v,RC, b,ib,jb)
		! (a,ia,ja) CSR formart sparse matrix
		! Rc='R' or 'C'. 'R' represents row; 'C' represents column
		real(kind=8), dimension(:), intent(in) :: a 
		integer, dimension(:), intent(in) :: ia, ja, v
		character, intent(in) :: RC
		real(kind=8), dimension(:), allocatable :: b
		integer(kind=8), dimension(:), allocatable :: ib,jb

		integer, dimension(:), allocatable :: vv, itmp
		integer :: ii, m, sv, jj, sib, kk, nnz

		m = size(ia) - 1
		sv = size(v)
		allocate(vv(sv))
		vv = v
		call qsort(vv,sv,kind(v),mso_cmp)	! sort vv
		! delete repetitive numbers
		if ( sv.gt.1 ) then
			ii = 2
			do while ( ii.le.sv .and. vv(ii-1).ne.vv(ii) )
				ii = ii + 1
			end do
			jj = ii + 1
			do while ( jj.le.sv )
				do while ( jj.le.sv .and. vv(jj-1).eq.vv(jj) )
					jj = jj + 1		
				end do
				if ( jj.le.sv ) then
					vv(ii) = vv(jj)
					ii = ii + 1
					jj = jj + 1
				end if
			end do
		end if
		sv = ii - 1

		if ( RC.eq.'R' .or. RC.eq.'r' ) then
			! delete row numbers>m
			ii = sv
			do while ( ii.ge.1 .and. vv(ii).gt.m )
				ii = ii - 1			
			end do
			sv = ii
			if ( sv.eq.0 ) then
				allocate(b(size(a)))
				b = a
				allocate(ib(size(ia)))
				ib = ia
				allocate(jb(size(ja)))
				jb = ja
				goto 100
			end if
			sib = m - sv + 1	! size(ib)
			allocate(ib(sib))
			allocate(itmp(sib-1))
			jj = 1
			kk = 1
			nnz = 0
			! value b, ib, jb
			do ii = 1, m, 1	
				if ( jj.le.sv .and. ii.eq.vv(jj) ) then
					jj = jj + 1
				else
					itmp(kk) = ii
					ib(kk) = nnz + 1
					nnz = nnz + ia(ii+1) - ia(ii)
					kk = kk + 1
				end if
			end do
			ib(kk) = nnz + 1
			allocate(jb(nnz))
			allocate(b(nnz))
			forall ( ii = 1:sib )
				jb(ib(ii):ib(ii+1)-1) = ja(ia(itmp(ii)):ia(itmp(ii)+1)-1)
				b(ib(ii):ib(ii+1)-1) = a(ia(itmp(ii)):ia(itmp(ii)+1)-1)
			end forall
		else
		end if 
	
		100 continue
		if ( allocated(itmp) ) deallocate(itmp) 
		if ( allocated(vv) ) deallocate(vv) 
	end subroutine mso_dcsrdltRC

	function mso_cmp(a,b)
		integer, intent(in) :: a,b
		integer(kind=2) :: mso_cmp
		if ( a.lt.b ) mso_cmp = -1
		if ( a.eq.b ) mso_cmp = 0
		if ( a.gt.b ) mso_cmp = 1  
	end function mso_cmp

	subroutine mso_prtspmtxz(a,ia,ja,fmt)
		implicit none
		complex(kind=8), dimension(:), intent(in) :: a
		integer, dimension(:), intent(in) :: ia, ja
		character(len=*), intent(in) :: fmt

		integer :: m, ii, jj

		m = size(ia) - 1
		do ii = 1, m, 1
			do jj = ia(ii), ia(ii+1)-1, 1
				write(*,'(I10,I10,10X'//'("("'//trim(fmt)//'","'//trim(fmt)//',"i)"))') ii, ja(jj), a(jj)
			end do
		end do
	end subroutine

	subroutine mso_prtSpMtxR(a, ia, ja, rows, cols, fmt, trans)
		use printVecMtx
		implicit none
		real(kind=8), dimension(:), intent(in) :: a
		integer, dimension(:), intent(in) :: ia, ja
		integer, intent(in) :: rows(2), cols(2)
		character(len=*), intent(in), optional :: fmt
		character(len=1), intent(in), optional :: trans

		integer :: i, mn(2), m, job(6), info, n
		real(kind=8), dimension(:,:), allocatable :: adns

		m = size(ia) - 1
		n = maxval(ja)
		allocate(adns(m,n))

		job = 1
		job(4) = 2
		call mkl_ddnscsr(job,m,n,adns,m,a,ja,ia,info)

		if ( present(trans) .and. trans .eq. 'T' ) then
			do i = cols(1), cols(2)
				call prtVecR( adns(rows(1):rows(2),i), fmt )
			end do
		else
			if ( present(fmt) ) then
				do i = rows(1), rows(2)
					call prtVecR( adns(i,cols(1):cols(2)), fmt )					
				end do
			else
				write(*,*) adns
			end if 
		end if
		if ( allocated(adns) ) deallocate(adns) 
	end subroutine mso_prtSpMtxR

	subroutine mso_dgmres_precondition(maxiter, relres, a,ia,ja,rhs, computed_solution, itercount)
		!*******************************************************************************
		!  Content:
		!  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
		!*******************************************************************************
		implicit none
		integer, intent(in) :: maxiter
		real(kind=8) :: relres
		real(kind=8), dimension(:), intent(in) :: a
		integer, dimension(:), intent(in) :: ia, ja
		real(kind=8), dimension(:), intent(in) :: rhs
		real(kind=8), dimension(:), intent(out) :: computed_solution
		integer, intent(in) :: itercount
		

    	integer :: n, m
    	integer :: ipar(128)
    	real(kind=8) :: dpar(128)
    	real(kind=8), dimension(:), allocatable :: tmp, b, residual
    	
    	integer :: rci_request, i
    	real(kind=8) :: dvar
    	real(kind=8), external :: dnrm2
    	real(kind=8), dimension(:), allocatable :: bilu0, trvec
    	integer :: ierr, err
    
    	n = size(ia) - 1	! the size of the problem

    	allocate( b(n))
    	allocate( residual(n))

		!---------------------------------------------------------------------------
		! Initialize variables and the right hand side through matrix-vector product
		!---------------------------------------------------------------------------

    	call dcopy(n, rhs, 1, b, 1)

    	computed_solution = 1.d0

    	ipar(15) = min(150,n,maxiter)
    	allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

    	call dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999
    	
    	ipar(8) = 1
    	ipar(15) = min(150,n,maxiter)
    	ipar(5) = ipar(15)
    	ipar(11) = 0
    	dpar(1) = relres

    	call dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
    	if (rci_request .ne. 0) go to 999

    	ipar(31) = 1
	    dpar(31) = 1.D-20
	    dpar(32) = 1.D-16
	    allocate(bilu0(size(a)), stat=err)
	    if (err /= 0) print *, "bilu0: Allocation request denied"
	    call dcsrilu0(n, a, ia, ja, bilu0, ipar, dpar, ierr)
		allocate(trvec(n))

	    if (ierr .ne. 0) then
	        write(*,'(a,a,i1)') ' error after calculation of the',' preconditioner dcsrilu0',ierr
	          go to 999
	    end if

    	1 call dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
    	
    	if (rci_request .eq.0) go to 3

    	if (rci_request .eq. 1) then
    		call mkl_dcsrgemv('n',n, a, ia, ja, tmp(ipar(22)),tmp(ipar(23)))
            go to 1
    	end if

    	if (rci_request .eq. 2) then
			! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
			! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, computed_solution, b, rci_request, ipar,dpar, tmp, itercount)
			! compute the current true residual via mkl (sparse) blas routines
            call mkl_dcsrgemv('n', n, a, ia, ja, b, residual)
            call daxpy(n, -1.0d0, rhs, 1, residual, 1)
            dvar = dnrm2(n, residual, 1)
            if (dvar .lt. relres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            call mkl_dcsrtrsv('l','n','u',n,bilu0,ia,ja,tmp(ipar(22)),trvec)
            call mkl_dcsrtrsv('u','n','n',n,bilu0,ia,ja,trvec,tmp(ipar(23)))
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
            !---------------------------------------------------------------------------
			! if rci_request=anything else, then dfgmres subroutine failed
			! to compute the solution vector: computed_solution(n)
			!---------------------------------------------------------------------------
        else
            go to 999
        end if

    	3 IPAR(13) = 0
    	call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)


    	if ( allocated(tmp) ) deallocate(tmp)
    	if ( allocated(bilu0) ) deallocate(bilu0)
    	if ( allocated(trvec) ) deallocate(trvec)
    	if ( allocated(b) ) deallocate(b)
    	if ( allocated(residual) ) deallocate(residual)  
    	return

    	999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
    	call mkl_free_buffers
    	stop 1
	end subroutine mso_dgmres_precondition

	subroutine mso_dpardiso(mtype,msglvl,a,ia,ja,b,x)
		USE mkl_pardiso
		IMPLICIT NONE
		integer, parameter :: dp = kind(1.d0)
		integer, intent(in) :: mtype, msglvl
		integer, dimension(:), intent(in) :: ia, ja
		real(kind=dp), dimension(:), intent(in) :: a
		real(kind=dp), dimension(:), intent(inout) :: b
		real(kind=dp), dimension(:), intent(out) :: x	

		!.. internal solver memory pointer 
		type(mkl_pardiso_handle), allocatable  :: pt(:)
		!.. all other variables
		integer :: maxfct, mnum, phase, n, nrhs, error, nnz
		integer :: error1
		integer, allocatable :: iparm( : )
		integer :: i, idum(1)
		real(kind=dp) :: ddum(1)
		!.. Fill all arrays containing matrix data.
		n = size(ia) - 1
		nnz = size(a)
		nrhs = 1
		maxfct = 1 
		mnum = 1
		!..
		!.. set up pardiso control parameter
		!..
		allocate(iparm(64))

		do i = 1, 64
		   iparm(i) = 0
		end do

		iparm(1) = 1 ! no solver default
		iparm(2) = 2 ! fill-in reordering from METIS
		iparm(4) = 0 ! no iterative-direct algorithm
		iparm(5) = 0 ! no user fill-in reducing permutation
		iparm(6) = 0 ! =0 solution on the first n components of x
		iparm(8) = 2 ! numbers of iterative refinement steps
		iparm(10) = 13 ! perturb the pivot elements with 1E-13
		iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
		iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
		iparm(14) = 0 ! Output: number of perturbed pivots
		iparm(18) = -1 ! Output: number of nonzeros in the factor LU
		iparm(19) = -1 ! output: mflops for lu factorization
		iparm(20) = 0 ! output: numbers of cg iterations

		error  = 0 ! initialize error flag
		! 		msglvl = 1 ! print statistical information
		! 		mtype  = -2 ! symmetric, indefinite

		!.. initialize the internal solver memory pointer. this is only
		! necessary for the first call of the pardiso solver.

		allocate (pt(64))
		do i = 1, 64
 		  pt(i)%dummy =  0 
		end do

		!.. back substitution and iterative refinement
		iparm(8) = 2 ! max numbers of iterative refinement steps
		phase = 13 ! 33 ! only solving

		call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, b, x, error)
		! write(*,*) 'solve completed ... '
		if (error /= 0) then
   			write(*,*) 'the following error was detected: ', error
   			goto 1000
		end if
      
		1000 continue
		!.. termination and release of memory
		phase = -1 ! release internal memory
		call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

		if (allocated(iparm))   deallocate(iparm)

		if (error1 /= 0) then
   			write(*,*) 'the following error on release stage was detected: ', error1
   			stop 1
		endif

		if (error /= 0) stop 1
	end subroutine mso_dpardiso




	subroutine mso_zpardiso(mtype,msglvl,a,ia,ja,b,x)
		USE mkl_pardiso
		IMPLICIT NONE
		integer, parameter :: dp = kind(1.d0)
		integer, intent(in) :: mtype, msglvl
		integer, dimension(:), intent(in) :: ia, ja
		complex(kind=dp), dimension(:), intent(in) :: a
		complex(kind=dp), dimension(:), intent(inout) :: b
		complex(kind=dp), dimension(:), intent(out) :: x	

		!.. internal solver memory pointer 
		type(mkl_pardiso_handle), allocatable  :: pt(:)
		!.. all other variables
		integer :: maxfct, mnum, phase, n, nrhs, error
		integer :: error1
		integer, allocatable :: iparm( : )
		integer :: i, idum(1)
		real(kind=dp) :: ddum(1)
		!.. Fill all arrays containing matrix data.
		n = size(ia) - 1
		nrhs = 1
		maxfct = 1 
		mnum = 1
		!..
		!.. set up pardiso control parameter
		!..
		allocate(iparm(64))

		do i = 1, 64
		   iparm(i) = 0
		end do

		iparm(1) = 1 ! no solver default
		iparm(2) = 2 ! fill-in reordering from METIS
		iparm(4) = 0 ! no iterative-direct algorithm
		iparm(5) = 0 ! no user fill-in reducing permutation
		iparm(6) = 0 ! =0 solution on the first n components of x
		iparm(8) = 2 ! numbers of iterative refinement steps
		iparm(10) = 13 ! perturb the pivot elements with 1E-13
		iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
		iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
		iparm(14) = 0 ! Output: number of perturbed pivots
		iparm(18) = -1 ! Output: number of nonzeros in the factor LU
		iparm(19) = -1 ! output: mflops for lu factorization
		iparm(20) = 0 ! output: numbers of cg iterations

		error  = 0 ! initialize error flag
		! 		msglvl = 1 ! print statistical information
		! 		mtype  = -2 ! symmetric, indefinite

		!.. initialize the internal solver memory pointer. this is only
		! necessary for the first call of the pardiso solver.

		allocate (pt(64))
		do i = 1, 64
 		  pt(i)%dummy =  0 
		end do

		!.. back substitution and iterative refinement
		iparm(8) = 2 ! max numbers of iterative refinement steps
		phase = 13 ! 33 ! only solving

		call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, b, x, error)
		! write(*,*) 'solve completed ... '
		if (error /= 0) then
   			write(*,*) 'the following error was detected: ', error
   			goto 1000
		end if
      
		1000 continue
		!.. termination and release of memory
		phase = -1 ! release internal memory
		call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

		if (allocated(iparm))   deallocate(iparm)
		if (allocated(pt))   deallocate(pt)

		if (error1 /= 0) then
   			write(*,*) 'the following error on release stage was detected: ', error1
   			stop 1
		endif

		if (error /= 0) stop 1

	end subroutine mso_zpardiso

end module MySparseOperator